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A COMPUTER PROGRAM TO CALCULATE ARTIFICIAL 
RADIATION BELT DECAY FACTORS 

by 

E. G. Stassinopoulos 
Goddard Space Flight Center 


SUMMARY 

The computer code presented here is a FORTRAN II, Monitor- 
Controlled program for the IBM 7090-7094, and has been developed 
for calculating artificial radiation -belt decay factors to evaluate 
the time -dependent change in the intensity of the Starfish electrons 
injected into the upper atmosphere on July 9, 1962. The program 
(called subroutine DECAY) has been written in subroutine form 
as part of an orbital flux code for calculating omnidirectional, 
vehicle-encountered fluxes. Execution time for the program 
varies, but even in the worst case it is still insignificantly small 
(a fraction of a second). By means of this code, anticipated fluxes 
at any epoch may be predicted from known radiation levels at a 
given time. 
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A COMPUTER PROGRAM TO CALCULATE ARTIFICIAL 
RADIATION BELT DECAY FACTORS 


by 

E. G. Stassinopoulos 
Goddard Space Flight Center 


INTRODUCTION 

A significant loss of Starfish electrons has occurred since their injection into the upper at- 
mosphere on July 9, 1962. It is assumed that this artificial radiation- belt population decreases 
exponentially with time (References 1 and 2), thus, the remaining particle intensity may be ade- 
quately approximated by a decay factor of the form, 




(1) 


where r is the mean lifetime of electrons with energies greater than 1.2 Mev (Reference 3) and At 
is the specified time difference (t 2 - tj . If the flux at time t 2 is then related to that of time by 
this decay factor, the expression, 




® (t i ) exp 



( 2 ) 


is obtained which may be used readily in a computer calculation. 



Measurements of mean lifetime, r, (by Van Allen and Lin with INJUN III and by Bostrom and 
Williams of APL with 1963-38c, References 1 and 2) have shown, that the decay constant is mainly 
a function of the magnetic shell parameter, L, and is almost independent of the field strength, B . 


It should be noted, that the 1963- 38c measurements were taken at a later epoch than the 
INJUN in and hence are consistently lower (Table 1 and Figure 1). 
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Table 1 


Subroutine-DECAY Mean Lifetime Measurements of Artificial 
Radiation- Belt Electrons. 


Measure- 
ment No. 

Magnetic- Shell 
Parameter, L 

Magnetic Field 
Strength, B 

Mean Life- 
time, T 1 1 

Mean Life- 
time, rt 

(earth radii) 

(gauss) 

(days) 

(days) 

1 

(1.15) 


(0) 

(0) 

2 

1.20 

0.185-0.205 

130 

120 

3 

1.23 

0.170-0.205 

190 

165 

4 

1.30 

0.100-0.230 

320 

235 

5 

1.40 

0.165-0.210 

500 

390 

6 

1.45* 


565* 

445* 

7 

1.50 

0.175-0.215 

590 

460 

8 

1.55* 


585* 

445* 

9 1 

1.60 

0.180-0.225 

(500) 

360 

10 

1.65* 


275* 

225* 

11 

(1.70) 


(100) 

(100) 

12 

1.725* 


50* 

65* 

13 

1.75* 


25* 

40* 

14 

(1.80) 


(0) 

(0) 

All values in parenthe 

ses are estimates not contained in original information. 



* Values indicated were 

obtained by fitting the 7-curve (Figure 1) to the source and to the estimates. 


tTj measurements made by Van Allen and Lin with INJUN III. 
tr 2 measurements made by Bostrom and Williams with 1963-38c. 




METHOD OF COMPUTATION 

The subroutine first calculates the time 
differenceAt - t 2 - t v where t l9 contained in- 
ternally, defines the epoch of the electron grid 
used* and t 2 , an input quantity, pertains to a 
given later date, for which flux information is 
being sought. 

If At < 0, the subroutine will bypass the rest 
of the program and, since no decay is involved 
in t^at case, it will set the decay factor equal to 
unity and return to the main routine. 

If At > 0 the subroutine continues, it then 
initiates computation of the decay factor by 



Figure 1 —Function r (L) fitted to measurements by Van 
Allen and Lin (tjJ Bostrom and Williams for L < 1.6 
and to estimates by Hess for L > 1.6. 


■"Presently the code is programmed for the E8 grid of November, 1962. 
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establishing once, in its first passage, the r-to-L dependence in the form of a look-up table where 
t(L) = 0 for L > 1.8 (see Table 1). Subsequently, it tests the L of the orbital position under con- 
sideration and if 1.15 < L < 1.8 it proceeds to obtain r by a linear interpolation from both sources 
available (i.e., Van Allen and Lin, and Bostrom and Williams) and then either computes or defines 
the proper decay factor, exp(- At/r). 

Finally, the subroutine sets up suitable control parameters in order to avoid, during its re- 
petitive use by the Orbital Flux Code, recomputation of any quantity that is both invariant with 
respect to a given trajectory and permanently established for that orbit during the first passage. 
In this manner the subroutine also passes over the once initialized look-up table. 


USAGE AND LIMITATIONS 

No actual limitations exist with respect to computation of a decay factor by the subroutine. It 
is programmed to accept any epoch t 2 including those that fall on, or before, its critical or ef- 
fective date . 

Epoch t 2 must be given in the form of year, month and day, and in that order. The subroutine 
will test each quantity separately, in the same order, and convert into days if necessary. 

There is no difficulty in substituting another effective date t j (for different flux data) into the 
code as long as appropriate reference relations are included. 

Similarly, the mean lifetimes may be changed easily without major modifications by simply 
inserting new tables of r to L values between statements 551 and 550 and by adjusting the limits 
of statements 550 and Dimension (if necessary). 

Due to the IBM 7094 floating point limita- 
tions, it is to be expected that the decay factors 
will be zero whenever the ratio At A becomes 
greater than that value, which will make the 
characteristic exceed the storage capacity of the 
register. When this happens, underflow occurs. 

Figure 2 shows the relationship of time differ- 
ence, At , to mean lifetime, r, where the shaded 
area indicates combinations causing such under- 
flow. This however does not affect the results 
in any way. Contributions of such small mag- 
nitudes can be neglected. A sample computa- 
tion with the r t decay constant for various con- 
stant flux ratios («/>/<£ 0 ) is shown in Figure 3. 


TIME DIFFERENCE, At (days) 



Figure 2— Underflow relationship of mean lifetime r and 
time difference A t. (Shaded area shows combinations 
causing decay factors to become zero.) 
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The constant r 1 is that measured by Van Allen 
and Lin with INJUN III. 

CONCLUSION 

Whenever the orbital flux code is employed 
and its output is utilized, caution must be ex- 
ercised in regard to the following shortcomings 
and restrictions: 

1. The subroutine applies to Starfish elec- 
trons only. It does not consider or in- 
clude natural electrons.* 

2. The mean lifetime (and hence the flux) 
has been set equal to zero for L > 1.8. 

3. The decay constants r x and r 2 were ob- 
tained from measurement of E > 1.2-Mev 
electrons. If used for all energy elec- 
trons, the resulting spectra may not be 
very accurate (values possibly inflated 
for low energy fluxes).* 



MAGNETIC SHELL PARAMETER, Hearth radii) 


Figure 3— Contours of constant flux ratius <J>/<t> 0 as func- 
tions of L and t for Starfish electrons. (Each curve traces 
a given fixed percentage by which an initial flux is re- 
quested to decay.) 


4. Since increasing solar activity in the years ahead will greatly affect the lifetimes of the 
Starfish electrons, the use of this routine for epochs later than 1966 is not advisable.* 
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Appendix A 

Subroutine 


Arguments and Parameters, Flow 
Chart and Listing 
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1. Arguments and Parameters 
A. Input: 

IN - Control parameter; governs the calculation of At; may have values 0, 1 or 2; at 

the start it is set equal to zero in main routine and remains unaffected by it 
thereafter; it is redefined in subroutine; it performs as follows: 

0: attempts to compute At , 

1: At is well defined, hence it bypasses section computing At and proceeds to 
obtain decay factors, 

2: At cannot be defined, hence it bypasses section computing At , sets decay fac- 
tors equal to one. 

DL - Magnetic shell parameter L associated with orbital position under consideration. 

YEAR - Year of date for which decay factor is to be calculated. 

AMONTH - Month of date for which decay factor is to be calculated. 


DAY - Day of date for which decay factor is to be calculated. 

NOT - Control parameter; governs the initialization of the r-L tables; may be -1, 0, +1; 

at the start it is set equal to zero in main routine and remains unaffected by it 
thereafter; it is redefined in subroutine; it performs as follows: 

-1: skips initialization of tables 

0, 1: initializes tables 


B. Internal: 
T1 

T2 

T3 

T 

TAL(I) 

TA1(I) 

TA2(I) 

KK 

TAUl 

TAU2 


- year(s) in days ''i 

- month(s) in days > These arguments establish the time difference At - ( 1 2 - t ,) 

- day(s) > 

- time difference At in days, i.e., the sum of Tl, T2, T3 

- L ; values of the r-L tables, where i l, • ■ • 14 

- r, values of the r-L tables, where i - l, 14 

- r 2 values of the t-L tables, where i r 1. •••14 

i 

- Control parameter within L-testing loop; defines special cases. 

- decay constant r x 

- decay constant r 2 
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C. Output: 


DFCTR1 - decay factor with respect to r, 

DFCTR2 - decay factor with respect to r 2 

Note: All arguments are in floating point notation except the parameters IN, NOT, KK which are 
fixed point quantities. 
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n n n n 


3. Listing 


SUBROUTINE DECAY! I N * DL , YEAR , AMONTH , DA Y , DPCTR 1 ,DFCTR2 * NOT ) 

THE MEAN LIFE TIME TAU1 IS OBTAINED FROM MEASUREMENTS MADE WITH 
INJUN III 

THE MEAN LIFE TIME TAU2 IS OBTAINED FROM MEASUREMENTS MADE WITH 
1963 38C 

DIMENSION TAL(20) ,TA1<20) *TA2(20) 

I F ( IN - 1 ) 315 * 320 *310 
315 T 1 = 0. 

T 2 = 0. 

T 3 = 0. 

T = 0. 

I F ( YEAR - 1962. >310*100*105 
IQO T 1 = -365. 

GO TO 110 

105 T 1 = 365 • * ( YEAR - 1963.) 

110 I F ( AMONTH - 1 1 . ) 2 04 * 200 * 202 
200 T 2 = 364. 

GO TO 210 
202 T 2 = 394. 

GO TO 210 

204 MONTH = AMONTH 

GO TO (90*91, 92*93*94, 95,96*97, 98*99) , MONTH 

90 T 2 = 60. 

GO TO 210 

91 T 2 = 91. 

GO TO 210 

92 T 2 = 119. 

GO TO 210 

93 T 2 = 150. 

GO TO 210 

94 T 2 = 180. 

GO TO 210 

95 T 2 = 211. 

GO TO 210 

96 T 2 = 241. 

GO TO 210 

97 T 2 = 272. 

GO TO 210 

98 T 2 = 303. 

GO TO 210 

99 T 2 = 333. 

210 I F ( DAY - 1 . >300*300*305 
300 T 3 = 1. 

GO TO 290 
305 T 3 = DAY 
290 T=Tl+T2+T3 

I F ( T - 484. )295*296*296 
296 T = T + 1. 

295 IF ( T ) 310*310*320 
310 DFCTR1 = 1. 

DFCTR2 = 1. 

IN = 2 
RETURN 

320 I F ( NOT )550*551 *551 
551 TaL<1) = 1.15 
T AL ( 2 ) = 1.20 
T AL ( 3 ) = 1.23 
T AL ( 4 ) - 1.30 
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T AL ( 5 ) = 1.40 
TAU6) = 1.45 
TAU7) = 1.50 
T AL < 8 ) = 1.55 
TAL(9> = 1.60 
TAL ( 10 ) = 1.65 
T AL (11)= 1.70 
T AL ( 1 2 ) = 1.725 
TAL ( 1 3 ) = 1.75 
TaU14)= 1.80 
TA1 ( 1 ) = 0. 

TAK2) = 130. 

T A 1 C 3 ) = 190. 

T Al ( 4 ) = 320. 

TA1(5\ = 500. 

T A 1 ( 6 ) = 565. 

T Al ( 7 ) = 590. 

TA1 (8) = 585. 

T Al ( 9 ) = 500. 

TAKIO)-- 275. 

TA1(11)= 100. 

T Al < 1 2 ) = 50. 

T A 1 I 1 3 ) = 25. 

TA1( 14)= 0. 

T A2 < 1 ) = 0. 

T A2 ( 2 ) = 120. 

T A2 ( 3 ) = 165. 

TA2(4)= 235. 

TA2(5)= 390. 

TA2(6)= 445. 

T A2 ( 7 ) = 460. 

TA2(8J- 445. 

T A2 ( 9 ) = 360. 

TA2< 10)=225. 

T A2 ( 1 1 ) = 1 00 • 

TA2< 12)= 65. 

TA2U3)= 40. 

T A2 ( 14 ) = 0. 

NOT = -1 

550 DO 900 K= 1 * 14 
KK = K 

IFITALIK) - DL ) 900 *500*501 
‘900 continue 
GO TO 502 

501 IFIKK-1 )5Q2 *502*321 

502 DFCTR1 = 0. 

DFCTR2 = 0. 

GO TO 370 

500 TAU1 = TAKK.M 
T AU2 = TA2<KK> 

GO TO 322 

321 TAU1 = TAl(KK)-( (TAl(KK)-TAl(KK-l) )/(TAL(KM-TAL(KK-in )*<TAL(KM 
1DL ) 

TAU2 = TA2t^M-t ITA21KK>-TA21KX-1> ) / * T AL l KK ) -T AL ( K.K- 1 ) ) ) * { T AL ( KX ) 
1DL ) 

322 DFCTR1 = EXPF ( - T/T AU 1 ) 

DFCTR2 = EXPFI-T/TAU2) 

370 I N= 1 

RETURN 

END 


Appendix B 


List of Symbols 


At = Time difference in days = (t 2 - 

t j = November 1, 1962, effective date relating to epoch of electron data 
t 2 = Any other date 

r = Mean lifetime of Starfish electrons in days 

Tj = Mean lifetime by Van Allen and Lin from measurements with INJUN III 
r 2 - Mean lifetime by Bostrom and Williams from measurements with 1963- 38c 
L = Magnetic shell parameter in earth radii 
B = Magnetic field strength in gauss 

<t> = Omnidirectional flux in particles per square centimeter per second (particles/cm 2 -second) 


NASA- Langley, 1965 G-629 
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